function state = rk4(func, statep, alphad, bankd, h)
    k1 = func(statep, alphad, bankd);
    k2 = func(statep + 0.5 * k1, alphad, bankd);
    k3 = func(statep + 0.5 * k2, alphad, bankd);
    k4 = func(statep + k3, alphad, bankd);
    state = statep + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
end